Dissipative Quantum Transport in Macromolecules: An Effective Field Theory 

Approach 
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We introduce an atomistic approach to the dissipative quantum dynamics of charged or neutral 
excitations propagating through macromolecular systems. Using the Feynman- Vernon path integral 
formalism, we analytically trace out from the density matrix the atomic coordinates and the heat 
bath degrees of freedom. This way we obtain an effective field theory which describes the real-time 
evolution of the quantum excitation and is fully consistent with the fluctuation-dissipation relation. 
The main advantage of the field-theoretic approach is that it allows to avoid using the Keldysh 
contour formulation. This simplification makes it straightforward to derive Feynman diagrams 
to analytically compute the effects of the interaction of the propagating quantum excitation with 
the heat bath and with the molecular atomic vibrations. For illustration purposes, we apply this 
formalism to investigate the loss of quantum coherence of holes propagating through a poly(3- 
alkylthiophene) polymer. 
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The investigation of the real-time dynamics of charged and neutral quantum excitations propagating through 
macromolecular systems is receiving a growing attention due to its potentially countless implications in nano-scale 
(opto-) electronics and in biophysics. For example, the study of the electric conduction through inorganicP^, organicP^ 
and biologicaP^L^ polymeric systems and aggregates is motivated by the perspective of realizing nano-scale, or 
even single-molecule^ transistors. In addition, the recent experimental observation of coherent exciton transport in 
photosynthetic protein-pigment complexes 13 has triggered a huge activity aimed at clarifying the interplay between 
quantum coherence, transfer efficiency and environment-driven noise^"^. 

Quantum transport processes in biomolecules and organic materials have been extensively discussed in the frame- 
work of phenomenological models in which the dynamics of the quantum excitation is described at the level of a 
one-body Hamiltonian, while the fluctuation-dissipation generated by the mole cular vibrations and by the heat bath 
are collectively represented by some effective bosonic fields (see e.g. Ref! n * 14 E^and references therein). The frequency 
spectrum of these bosons is modeled phenomenologically, encoding information gained from MD simulations and ex- 
O periment. These effective models provide computationally efficient tools to study the global and general mechanisms 
which underlie the long-range charge transport in macromolecules and investigate de-coherence and re-coherence 
^^ phenomena. 

^ In order to establish a more direct connection between the quantum transport dynam ics and the specific physico- 

chemical properties of the molecule under consideration, an alternative approac h 19 * 1 ^^ has been developed in which 
the atomic coordinates are treated explicitly and are evolved in time using a classical MD algorithm. The parameters 
of the effective tight-binding Hamiltonian are derived from the electronic structure, hence depend on the instantaneous 
molecular configuration. They are evaluated at periodic time intervals, within the Born-Oppenheimer approximation, 
e.g in the DFT-TB scheme, while the current flowing through the system is estimated in Landauer theory. This 
method neglects the effect of the coupling between the quantum excitation and the atomic nuclei on the molecular 
dynamics. In addition, the charge is assumed to propagate instantaneously across the system. In general, the validity 
of these assumptions may be questionable whenever the quantum transfer dynamics and the molecular dynamics 
occur at comparable time scales. 

A theoretical framework which takes into account of all the correlations between the quantum excitation dynamics 
and the molecular dynamics is the so-called non-equilibrium Green's function (NEGF) method (see e.g. RefP^ and 
references therein). This formalism was introduced to describe electric conduction across molecular wires and has 
been extensively applied to systems consisting of a few hundreds of atoms. In the NEGF approach, the interaction 
with the metallic leads and with phonons is treated at the quantum level, since the method is often used to investigate 
conduction at very low temperatures. The real-time dynamics of the system is described in the Keldysh formalism by 
a Dyson equation containing different types of Green's functions, corresponding to different segments of the Keldysh 
contour. Unfortunately, this feature makes the approach quite cumbersome and computationally demanding, hence 
hardly applicable to large molecular systems, such as DNA or of conjugate polymers. 

Fortunately, when investigating the transport of electric charge or neutral quantum excitations across macro- 
molecules at room temperature, a quantum description of the molecular vibrations is not really mandatory, as it is 
demonstrated by the success of classical molecular dynamics (MD) simulations. On the other hand, in these studies 
it is important to take into account of the fluctuation and dissipation generated by the solvent (in biomolecules) or 
by neighboring molecules ( in organic frameworks ) . The NEGF approach does not exploit the possibility of taking 



the classical limit for the molecular vibrations and does not explicitly account for the fluctuation-dissipation effects 
induced by the heat bath. 

In the present paper, we develop a formalism to describe quantum transport in macromolecular systems at room 
temperature, in which we exploit the possibility of treating the molecular vibrations at the classical level. In analogy 
with the NEGF method, the equations of motions for the entire system are not postulated phenomenologically, but 
rather rigorously derived from the reduced density matrix of the system. This way, one retains in a consistent way 
the relevant couplings between the molecular, heat bath and quantum charge degrees of freedom. 

The computational efficiency of our approach is enhanced by the fact that we are able to analytically integrate 
out the vibronic and heat bath dynamics and obtain a simpler (but rigorous) effective theory, which is formulated 
solely in terms of the quantum excitation degrees of freedom. In this theory, fluctuation-dissipation effects and 
thermal oscillations of the molecule are taken into account through effective interaction terms, which are derived from 
first principles. The information about the configuration-dependent electronic structure of the molecule is implicitly 
encoded in the parameters which appear in the effective theory. These are determined once and for all by means of 
quantum-chemistry calculations. 

A second important feature of our theory is that it allows to strongly simplify the formalism to describe the non- 
equilibrium dissipative dynamics. Indeed, our final path integral representation of the density matrix is formally 
equivalent to that of a vacuum-to- vacuum Green's function in a zero-temperature quantum-field theory. In particular, 
the time variable in the action functionals is integrated along the real axis, and not along the Keldysh contour, like 
in the NEGF approach. The formal connection with quantum field theory in vacuum makes it straightforward to 
derive Feynman diagrams to perturbatively compute the effects of the dissipative coupling between the propagating 
quantum excitation, the heat bath and atomic degrees of freedom. These corrections are obtained analytically, i.e. 
without any need to average over many independent MD trajectories. 

As a first illustrative application of this formalism, we develop a simple model to describe intra-chain charge prop- 
agation in a poly(3-alkyothiophene) (P3HT). First, all the couplings in the effective theory are derived from quantum 
chemistry calculations, at the DFT-TB level. Next, the leading-order perturbative corrections to the evolution of the 
charge density distribution are obtained by computing a few Feynman diagrams. The results are then compared with 
those obtained by solving numerically the coupled quantum/stochastic equations of motion for the polymer configu- 
rations and the charge reduced density matrix. From this comparison, we are able to assess the range of parameters 
and time intervals over which we expect the perturbative approach to be reliable. We also monitor the progressive loss 
of quantum coherence and identify the effective interactions which are responsible for de-coherence and re-coherence 
phenomena. 

The manuscript is organized as follows: The formalism and path integral representation of the density matrix is 



developed in sections O and III In section IV we derive the perturbation theory and compute the relevant-leading 



order Feynman diagrams. In section M we present our illustrative application to the P3HT system. Conclusions and 



perspective developments are summarized in section VI 



II. MODELING THE DYNAMICS OF QUANTUM EXCITATIONS IN MACROMOLECULES 

An entirely ab-initio approach to quantum transport processes in macromolecules would involve solving the complete 
time-dependent Schodinger equation which couples all the nuclear and electronic degrees of freedom, both in the 
molecule and in its surrounding environment. Clearly, such an approach is beyond the reach of any present or 
foreseen computational scheme. 

A commonly adopted framework to reduce the computational complexity of this problem consists in relying on the 
Born-Oppenheimer approximation and in coarse-graining of the electronic problem into an effective configuration- 
dependent one-body tight-binding Hamiltonian. In such an approach, the molecule is first partitioned into several 
fragments, hereby labeled by an index m, each of which is assigned a frontier orbital |0 m ). Such a partition of the 
molecule must be defined in such a way that the electron density is significantly more delocalized within each molecular 
fragment than over different fragments^. The frontier orbitals \(/> m ) are calculated by solving the Schrodinger equation 
for a reduced portion of the molecule, centered at the fragment m. The system's wave function is then obtained by 
diagonalizing the Hamiltonian projected onto the space of the frontier orbitals |0 m ). 

For example, in studying electronic hole transport in DNA, the molecular fragments can be chosen to coincide with 
the base-pairs and their highest occupied molecular orbitals (HOMO's) can be obtained by solving the Schrodinger 
equation for an isolated paii^M The propagation of the charge carriers through the molecule is then modeled as the 
hopping of holes between neighboring base-pairs. 

In general, it is convenient to choose the dimensionality of the molecular fragment index m according to the topology 
of the system under consideration. For example, in molecular wires it is natural to adopt a one-dimensional index (see 
left panel in Fig. HI). On the other hand, if the electron can also propagate along the side-chains of branched polymers, 
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FIG. 1: Labeling of molecular fragments in molecules with different topology. In linear polymers, the fragments among which 
the charge hops form an effective one-dimensional system ( left panel). In branched polymers, a two dimensional representation 
may be introduced in order to account for propagation along the side chains (center panel). In polymer assemblies a third 
index may be used to label the different molecules (right panel). 



a two-dimensional index is most appropriate (see center panel in Fig. pi). In assemblies of branched polymers, a third 
component of the index m may be introduced to distinguish between the different monomeric units. 

In the following, we shall denote with r % = (r£,r*,r*) the Cartesian coordinates of the i— th atom. The set of all 
atomic nuclear coordinates is collectively represented by the configuration space vector Q, 



The dynamics of the entire system is modeled by the following quantum Hamiltonian: 

H = Hmc + Hm + Hb + Hmb- 
In this equation, Hmc the tight-binding Hamiltonian for the quantum excitation, 

HmC — y , ) y /mn[0] ^L,s^n,s- 
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We stress that the transfer matrix elements depend on the molecular configuration Q. The a^ (a m?s ) operators 
create (annihilate) a quantum excitation at the molecular fragment m and s denotes the third component of the 
spin. For sake of definiteness, in this work we shall focus on the propagation of electron holes in Highest Occupied 
Molecular Orbitals (HOMO's), which is relevant for many molecular charge transfer phenomena. Correspondingly, 
the creation and annihilation operators are assumed to obey the anti-commutation relations 
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The generalization of the present formalism to the case in which the propagating excitation is bosonic (e.g. an 
exciton) is straightforward and will not be discussed explicitly. Depending on the specific molecular system, it may be 
necessary to include also the coupling between different molecular orbitals. This can be done by introducing additional 
creation- annihilation operators. 

Sometimes it is convenient to split the f mn [Q] matrix elements into the hopping matrix elements T mn [Q] and on-site 
energies e m [Q] S mn : 



/mn(Q) = T mn (Q)(l - 6 mn ) - e n (Q)5 n 



(5) 



The parameters T mn and e m are obtained from the fragment orbitals |0 m ) and |0 n ) and depend parametrically on 
the configuration vector Q: 



r m n(Q) 

e m (Q) 



v\Hel.\(/>m), 



(6) 
(7) 



where H e i is the electronic Hamiltonian (for example, the Kohn-Sham Hamiltonian of density functional theory). In 
a molecular wire, the Hamiltonian Hmc must include also the coupling with the leads which play the role of electron 
source and sink. 



The Hamiltonian Hm in Eq. (I2| governs the conformational dynamics of the atomic nuclei 1 , in the absence of 
electronic holes and reads 

H M =^ + V(Q), (8) 

where the P is the momentum canonically conjugated to the configuration vector Q. V(Q) is the molecular potential 
energy, evaluated in the Born-Oppenheimer approximation. This includes the interaction between the different atoms 
within the molecule and possibly with the external fields. We stress that the potential energy V(Q) in Eq. ^ 
depends only on the molecular configuration. This is equivalent to taking the adiabatic limit for the dynamics, and 
to assume that the location of the quantum excitation does not alter in a significant way the interaction between the 
atomic nuclei. In many cases of interest, the validity of this approximation has been questioned and corrections to 
the adiabatic regime have been proposed. However, in this first paper we shall not deal with these complications. 

The part of the Hamiltonian Hb + Hmb describes the coupling of the molecule with a thermal heat bath in the 
Leggett-Caldeira model 22 , i.e. through an infinite set of harmonic-oscillators coupled to each atomic coordinate: 

fijUjpU , (9) 

3N oo / 2 \ 




Hmb = J2J2[ -CjXjQ* + %^J2& • ( 10 ) 

OL=l j = l \ ^ 3 / 

X = (xi,^2, . . .) and n = (7Ti,7T2, . . .) are the harmonic oscillator coordinates and momenta, fij and uoj denote their 



masses and frequencies and Cj are the couplings between atomic and heat bath variables. The last term in Eq. (10) 
is a standard counter-term introduced to compensate the renormalization of the molecular potential energy which 
occurs when the heat bath variables are traced out (see e.g. the discussion in RefPS). 

The model introduced so far represents the starting point of many approaches which have been proposed to describe 
quantum transport in molecular systems. For example, in the method used in RefPand Ref. 21 , the evolution of the 
molecular degrees of freedom Q is described at the classical level by means of a MD algorithm based on the force fields 
obtained by parametrizing the molecular potential energy V(Q). Then, the current flowing through the molecule is 
evaluated at regular time intervals, in the Landauer formalism, using the instantaneous values of the f\ m [Q(t)] tight- 
binding coefficients. Such an approach retains the effects of the molecular thermal oscillations on the charge dynamics. 
Indeed, fluctuations on Q generate dynamical disorder in the tight-binding matrix elements f\ m (Q). On the other 
hand, this approach neglects the back-reaction of the charge dynamics on the molecular dynamics. For example, 
it does not account for the fact that the system's total energy may be decreased by visiting some specific atomic 
configurations in which some on-site energies of the quantum excitation are lower. 

In a recent paper 2 ^, the path integral formalism was used to explicitly eliminate the dynamics of the heat bath 
variables, and take the classical limit for the atomic degrees of freedom. As a result, an algorithm was derived to 
describe in a dynamical way the coupled evolution of the molecule and the charge in the heat bath. In particular, 
in the high-friction limit, the quantum and stochastic evolution of the charge wave function |\£) and the molecular 
configuration Q in an infinitesimal time interval At is determined by the following equations: 



q a (t + At) = q a (t) - ^^ (V[Q(t)\ + £ m > m »(*) / nm (Q)]) + J2™^) 

(11) 

|tf(t + At)) =e~^ HMc[Q(t)]\^(t)), 

where p\ m (t) = (\I/(£)|l)(m|\I/(£)) is the time-dependent (reduced) density matrix of the hole, 7 is the heat bath 
friction coefficient and £ a (t) is a stochastic variable sampled from a Gaussian distribution with zero average and 
unitary variance. Since the motion of the molecular degrees of freedom is stochastic, the charge probability density 
at different times is obtained from the average over many independent trajectories, which may turn out to be a 
computationally challenging procedure. 

In the next section, we review and further develop the path integral approach to hole-transport in macromolecules 
and obtain a scheme which does not require to perform any MD or Langevin simulation. 



1 Notice that, for sake of notational simplicity, we are assuming that all atoms have the same mass M. The generalization to different 
atomic masses is straightforward and will not be discussed here. 



III. EFFECTIVE FIELD THEORY FOR THE REDUCED DENSITY MATRIX 



Let us assume that a hole is initially created at the HOMO of some molecular fragment Iq. We are interested in 
computing the conditional probability Pt(k/, |lq) for the hole to be found at the HOMO of the molecular fragment 
kj, after a time interval t. Such a probability is described by the following time-dependent reduced density matrix 
element: 



^(k/IkO 



TV[ Ik/Xk/lflt) ] _ TrJlkfHkfle-^ p(0) e*™] 



Tr p{t) Tr ,3(0) 

where p(0) is the initial density matrix, which is taken to be in the factorized form 



/3(0) = IkiXki 
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(12) 



(13) 



Eq. ( 13) corresponds to assuming that, at the initial time, the molecular degrees of freedom and the heat bath degrees 



of freedom can be considered separately equilibrated at the same temperature. Hence, the normalization factor at 
the denominator reads 



IV p(0) = 2 • Z M (P) xZ b (/3), 



(14) 



where Zm and Zb are the quantum partition functions for the molecule and heat bath degrees of freedom and the 
degeneracy factor 2 follows from enumerating the initial spin states. We emphasize that the ratio in Eq. (12) expresses 
dynamics of all the degrees of freedom at the fully quantum level. 

Let us now derive a path integral expression for its numerator. To this end, we choose to adopt a second-quantized 
representation 2 of the quantum charge dynamics, while retaining the standard first- quantized representation for the 
dynamics of the atomic coordinates and for the harmonic oscillators in the heat bath. 

The path integral representation of the reduced density matrix (12) is obtained by performing the Trotter decom- 



position of the forward- and backward- time evolution operators e %Ht and e lHt and of the imaginary-time evolution 



operators e k b 



-Ha 



and e 



-Hb 



which appear in Eq.s (12) and (13). In practice, our choice of the representation 



of the charge, the heat bath and the molecule dynamics corresponds to introducing the following resolution of the 
identity: 
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(15) 
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where the states |Q,X, 3>) collect the set of all hole's coherent states (constructed from the annihilation operators 
associated to each molecular fragment, ak, s ), the set of the eigenstates of the molecular coordinate operator Q and 
the set of eigenstates of the heat bath generalized coordinate operator X. Throughout this paper we shall adopt 
Einstein's notation and implicitly assume the summation over all repeated indexes, except for the initial and final 
position of the molecule, Iq and kj, which are held fixed. 

Once the conditional probability (12) is written in the path integral form and the Gaussian functional integral over 
the harmonic oscillator variables is carried out analytically, one reaches the expression 
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(16) 



2 As we shall see below, the field-theoretic description for the quantum charge is adopted in order to obtain a simpler final representation 
of the density matrix. Indeed, it allows to replace the time integration along the Keldysh contour (shown in Fig. [2} with a standard 
time integral along the real axis. 



L 3m t 
A 








Q'(0) = Qi 




5R( 






V 








I 






Q(0) = Q"(0) = Q 


Q'(t)=Q"{t) = Q f 




-03 


Q(-iP) = Q % 







FIG. 2: Boundary conditions of the molecular configurations paths defined on the Keldysh contour appearing in the path 
integral (16). 



where f dQf, J dQi denote the standard (i.e. Riemann) integral over the final and initial molecular configurations, 
respectively, J dQ is an integral over some intermediate configuration and VQ denotes the functional integral measure 
(see Fig. |2]). The functionals which appear at the exponent are defined as 
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B(t) is a Green's function which encodes the fluctuation-dissipation induced by the heat bath and reads: 
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The time scales at which thermal oscillations are damped and memory effects in the heat bath are relevant can be 
tuned by varying the frequencies of the virtual harmonic oscillators in Eq. (10). In particular, here we consider the 
so-called Ohmic bath limit, in which the B(t) reduces to 



B(t)^B ohm (t) 
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It is well known that, in the classical limit for the molecular motion, this choice for B(t) leads to a Langevin dynamics 
with friction coefficient 7 (see e.g. Ref. 23 ). Hence, Eq. (16) represents a quantum generalization of the stochastic 
dynamics of the molecule, which includes also the time-evolution of the hole. As usual, the path integral (16) is 



defined over Grassmann fields or complex fields, depending if the propagating excitation is a fermion or a boson. 

It is important to note that the Q (t) and Q (t) variables which appear in the path integral (16) represent the 
configuration of the molecule propagating forward and backwards in time respectively, while the Q(t) variables are 
associated to the evolution of the same molecule along the imaginary time direction. All these paths can be collectively 
represented by path variable integrated along the so-called Keldysh contour (see Fig. [2| . 

Equivalently, the path integral (16) can be expressed in a form in which forward and backward molecular paths 



Q'(t) and Q"(t) are replaced by their average and difference, respectively: 
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The result is 
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where the functionals at the exponent are defined as 
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In the path integral representation (16), the time-evolution of the charge- molecule system in contact with a dissipa- 
tive heat bath follows directly from the quantum Hamiltonian defined in Eq. pL without any further approximation. 
In particular, the molecule's configurational dynamics and the charge's quantum propagation are described at the 
fully quantum level. In addition, there is no restriction either on the strength of the coupling between the molecu- 
lar vibrations and the charge, nor on the amplitude of the conformational changes which the molecule can undergo 
within the time interval t. Clearly, computing such a path integral is a formidable task and further approximations 
are needed. 

Our first approximation consists in taking the classical limit for the dynamics of the molecular atomic coordinates. 
To this end, we begin by noting that the saddle-point equations which are derived by functionally differentiating the 
exponent in Eq. (16) with respect to R,y,(j) > \(j>" lead to the condition y(t) = for all t (see derivation in Ref!^). 
Following the discussion in RefP^, we impose the classical limit on the molecular motion by assuming that the 
path y(t) remains close to its saddle-point configuration (hence represents a small fluctuation) and by imposing the 
boundary-condition y(0) = 0. 

We can verify that the correspondence principle is fulfilled by such an approximation. Indeed, in the absence of the 
quantum charge, the resulting expression for the conditional probability coincides with the Onsager-Machlup path 
integral representation of the Langevin dynamics of the molecule in its heat bath. To prove this, let us provisorily 
drop all the coherent fields, expand the functionals in the exponents to quadratic order in y and perform the resulting 
Gaussian integration. This way, the path integral reduces to 
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where Som[R] is the well-known Onsager-Machlup functional, which assigns a statistical weight to the stochastic 
trajectories in the Langevin dynamics: 
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Finally, we take the saddle-point approximation for the path integral DQ, which corresponds to taking the classical 
limit also for the partition function of the initial molecular configuration. The result is 
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where po(Qi) 
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is the initial distribution of molecular configurations. We recognize that this is the nor- 



malization condition on the solution of the Fokker-Planck equation, expressed in path integral for m 125 * 26 -! anc [ i s the 
starting point of the so-called dominant-reaction pathway approach to investigate the long-time dynamics of macro- 
molecules 2 ^^. 



Let us now return to the path integral expression in Eq. (23), in which a quantum charge is allowed to propagate 
across the molecule and discuss our second approximation. We note that the quantum transport dynamics is in 
general much faster than the characteristic time scales for major conformational transitions of macromolecular systems 
(typically ranging from few nanoseconds to many seconds or even larger). Hence, during the time intervals which are 
relevant for quantum propagation phenomena, the molecule can be assumed to follow at most only small oscillations 
around the mechanical equilibrium configuration Qo, which is defined as the global minimum of the molecular potential 
energy V(Q). In this small-oscillation limit, it is convenient to introduce the atomic displacement variables 

Sr(t) = R(t) - Q , (29) 

and regard both the Sr(t) and y(t) as small quantities of the same order. 

In the expansion of the W functional up to quadratic order in Sr and y we obtain a term 

= Snyjl-Lij + . . . , (30) 

where Hij = QQ d QQ. V(Q)\q=Qo is the Hessian matrix of the potential energy at the point of mechanical equilibrium. 
A small deviation from the equilibrium configuration Qo generates a small change in the hopping matrix elements 
and in the on-site energies which define the tight-binding Hamiltonian pi). To the leading-order in the Taylor expansion 
in powers of 5r and y we have: 

/nm(r-f) = /° m + /im(^i-|)+..., (31) 

/nm(r+|) = /° m + /Am(Wf )+.-., (32) 

where /„ m = / nm (Qo) an d /nm — ^ifnm(Q)\Q=Q - In the small oscillation regime, the path integral over y is of 
Gaussian type and can be performed analytically, yielding 
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where the functional integral over Sr(t) is unrestricted also at time and time t and the S functional reads 
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It is convenient to introduce a differential operator L, which depends on molecular coordinate indexes z, j: 

[L] y = M (d 2 t , + 7 ft/) Sij + Hij , (35) 

its Hermitian conjugate reads 

[L% = M (d 2 t , - 7 ft') S tj + Hij, (36) 

Note that in L and L^ the time- derivatives are defined to act on the right. Using such operators, the functional in 



Eq. (34) can be rewritten as 
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The path integral ( 33 ) describes the coupled dynamics of the nuclear coordinates and of the electronic hole in the 



molecule. It is instructive to first consider this conditional probability in the limit in which the couplings between the 
charge and the molecular degrees of freedom are completely neglected. In this case, the path integral factorizes as 



Ptfrfto 



In this equation, 



|G£(k/,*|ki)| 2 x 
|G£(k/,t|kO| 2 . 



1 



Z M (P) 



T)5r e iM ~< ■'° 



jk /o dT(M5r-~fSr+m 3 Sr 3 ) 2 - §Sr i (0)n ij Sr i (0) 



e 2 



G>(k f ,t\ki) = (k f \e-^ ot I ki,), with H = / lm (Q ) a\ a m , 



(38) 
(39) 



(40) 



is the hole's Green's function defined by the tight-binding Hamiltonian Hq •> which is evaluated keeping the molecule 
"frozen" in its minimum-energy configuration Qo- 



A. Dirac-Like Notation 



Let us now return to the most general case, in which the interactions between the hole, the molecule and the heat 
bath are fully taken into account. The symmetric structure of the exponent in the path integral representation of 
the conditional probability P(kf,£|k^) suggests to collect all coherent field degrees of freedom </>>,, </>,, </k, 0. into 
a single 4-compenent Grassmann field ip defined as: 



^n 



( ^.t \ 



(41) 



Similarly, we collect all conjugate fields into ip^ = ((/>* ^, (f)* , , 0* *., 0* . 1 . In view of the formal analogy with Dirac 

theory 3 it is convenient to introduce also the following 4x4 matrixes, which define the projection onto the upper and 
lower spinor components and the interchange between them: 



7o 



10 
10 
0-10 
0-1 



7+ 



/ 1 

10 



\o 0. 



/o 

10 

Vo 1 



75 



/o 1 
1 
10 

Vo 1 



(42) 



3 Clearly, in the case in which the propagating particle is a scalar boson, the ip field has only scalar field two-components, ip = ((f) 1 \<j)") 
and the gamma matrixes are replaced by Pauli matrices 70 — > T3 , 75 — > n . 
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In addition, we change variable in the integration over the tp^ field by means of the substitution 

Mt) = iPi(t) 70. (43) 



Using the notation defined in Eq.s (41), (42) and (43), the conditional probability is written as: 

x expQs [^,^]J exp(-S e// .[<Jr]) exp Q {h[5r^^] + I 2 [Sr,^] + J 3 [<M} J , (44) 



where the terms 



5 [^,^] = / dt' ^m (^^/ - /^ n ) ^n , 
*/0 

S eff [Sr] = ^-J**>Sr*[Li-L] 



(45) 



,5H' , (46) 

describe the quantum propagation of the charge and the Langevin dynamics the molecular coordinates, in the absence 
of any coupling between them. The functionals I\ , I 2 and I3 are the interaction terms and are defined as 

h[8r, $,%!)] = - J dt' f 5r l (47) 

Jo 

*<M = ^ f& 4 \L\ij W, (48) 

l3 ^^ ] = W^I* dt ' Jl ° Jl °' (49) 

with J % = i/j m f^ nn ijj n and Jq = ijj m 70 f^ nn ip n - We note that the couplings I\ and I 2 vanish in the classical limit 
■jjr- — >> 0. The surface terms C±(t,0) and £2(^0) follow from the over completeness of the coherent-field basis and 
from the Boltzman distribution of the initial configuration, respectively, and read 

AM) = (^ m (0) 70 7 + ^m(0) +^ m (t) 70 7- fcW) , (50) 

C 2 (t,0) = ^(0)^^(0) + ... (51) 



Some comments on Eq. (f44j) are in order. Firstly, we note that the overall minus sign appearing in front the 
integral is a consequence of the Fermi statistics and ensures the overall positivity of the probability density. Secondly, 
we observe that, while the path integral ( [33]) is defined over forward- and backward- propagating fields (i.e. along 
the Keldish contour), the path integral Eq. (|44|) contains only the integration in the forward time direction. Indeed, 
the backward-propagating fields have been replaced by lower-components of the "spinor" field, hence can be formally 
interpreted as anti-matter degrees of freedom propagating forward in time. 

A major simplification which follows from our approximations is that the integral over the small displacement of 
the molecular coordinates from their equilibrium position Sr can be performed analytically. The result is an effective 
theory with non-instantaneous interactions between the holes: 

Pt(k f \ki) = - J T>j, Vi/> e- Cl W (^ k/ (t) 7 _75^k / (t)V; k ,(0)7 + 75^(0)) e~^ s ^^ 

x e mfo dt ' dt " Jo(t')Vi 3 {t'-t")r{t")-f£ fi dt'dt"j i (t')A ij (t'-t")j i (t") , 52 x 

In this equation, &%j{t f — t"), V%j{t — t') are respectively the Green's functions of the [LM] operator and the sum of 
the Green's functions of the L and 77 operators. 

In order to explicitly compute them, it is convenient to transform to the normal mode basis, in which the Hessian 
of the potential energy at the minimum energy configuration Q is diagonal: 

ul n SJ u JZ = s kz m n 2 k . (53) 

In this equation, Q^ denotes the frequency of the fc-th normal mode. 
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In this basis, the expressions of the vibron propagators A^(£' — t") and Vij(t — t') read (see the derivation in the 



appendix \M 






e -h 1*1 

2 m 2 n 2 k 

2 e~h W 

M uj$ 



u lk U kj a>k(i), 



U} k U kj b k {t), 



(54) 
(55) 



where ujq = \/|4fi| — 7 2 | and 



a fe (i) 



b k (t) 



sin (I Wo fe |i|) 



cos (|wg |t|) 



sinh (|w§ |i|) cosh (I w g |t|) 



w 6 
sin (Jo,* |i|) 
sinh (| w g |t|) 



7 



if 20 fe > 7 , 
if 2f2 fe < 7 , 



if 2il k > 7 , 
if 2f2 fc < 7 . 



(56) 

(57) 
(58) 



It is useful to consider the asymptotic expressions for the Green's functions in the limit 7 3> 2il k , which corresponds 
to the over-damping regime: 



Ay(t) = 

Vij{t) = 



1*1 



2 M 2 7 fig 
1 



*4 ^ 



M7 



2^1 

7 



^1*1 



-i\t\ e -^\t\ 



ui u ki . 



In the opposite under-damped regime ( 7 <C il k ) t ne asymptotic expression for the propagators is: 



Ay(t) 
Vy(«) 



2 M 2 7 ft 2 
e -l 1*1 



cos(n fe |t|) [// fc t/fej, 



2 M il k 



sm(il k \t\) U} k U kj . 



(59) 
(60) 

(61) 
(62) 



Eq. (52) is one of the central results of this work. It shows that the conditional probability P(ky,£|lq) for the 
dissipative dynamics of the quantum charge can be written in a form which is formally analog to that of a vacuum-to- 
vacuum amplitude, in an effective zero-temperature Dirac-like quantum field theory. This analogy is quite remarkable, 
since our theory describes an open system and is fully consistent with the fluctuation-dissipation relation. We also 
emphasize that the path integral expression (44) is real and positive definitive as it yields directly the conditional 



probability, even though it is formally equivalent to probability amplitude (namely a two-point Green's function) in 
the effective quantum field theory. 

We remark again that a particularly attractive feature of this formulation is that it does not involve the Keldysh 
contour. This strongly simplifies the formalism, since one does not need to introduce different types of Green's 
functions to distinguish between different sectors of the Keldysh contour. Instead, one can adopt time-ordered 
Feynman propagators and apply text-book perturbative and non-perturbative quantum field theory techniques to 
evaluate directly the conditional probability and the expectation values of operators. 



IV. PERTURBATION THEORY AND FEYNMAN DIAGRAMS 



In the short-time and weak-coupling regimes, the conditional probability P t (k/|k^) can be computed analytically 
in a perturbation theory derived by performing a Taylor expansion of the exponents in Eq. ( 52 ) in powers of the 
interaction terms 



Vl = -S / dt ' l dt "^(t') f ma V-n(t') &ij(t'-t") 4> m ,(t") P mn Vw(f'), 

ph Jo Jo 



v > - s I? 1 !?" 



^m(i') 70 fLn V>n(0 Vy^ - t") ^ m ,(t") fU Vv (O • 



(63) 

(64) 
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(65) 



FIG. 3: The diagram corresponding to the unperturbed contribution to the conditional probability. 
The conditional probability is then written as 

oo 

p t (kflk i ) = ^p t (i) (k f lk i ), 

% 

where P\ (kf |ki) corresponds to the unperturbed conditional probability, which neglects all the couplings between 
the hole, the heat bath and the vibronic modes, 

P t (0) (k, |k,) = -^ J-Dj ZV e~ c ^ (0~ k/ (i) 7-75 lM*) ^(0) 7 + 75 lM0)) e~* s ^ . (66) 

Its normalization factor Z^ can be written in path integral form as: 



Z (0) 



D0 Vxj) e- £l (*'°) J2 OM*) 7-75 0k/ (*) 0k, (0) 7 + 75 0k, (0)) e~* 



iSo[^,^] 



(67) 



The leading-order perturbative correction in the series ( 65 ) reads 
-1 



P{ 1] (kf |kt) 



D0 2?V e- £l (*'°) (0 k/ (i) 7_ 75 k/ (t) kl (O) 7+75 0k,(O)) (Vi + K 2 ) e -* s °^, (68) 



Z(°) + ZC 1 ) 
where the corresponding leading-order correction to the normalization factor is 

Z« = \v^V^e- c ^m ^(0 k/ (i)7_ 75 ^ k/ (t)0 ki (O)7 + 75 0k,(O)) (Vi + V 2 ) e -* 



Sol*,*] 



(69) 



Eq.s (66), (67), (68) and (69) correspond to correlation functions in the free limit for the effective Dirac-like 



quantum field theory. According to Wick's theorem, these Green's functions can be evaluated by considering the sum 
of all possible contractions between the ip and ip fields and replacing each contraction with time-ordered Feynman 
propagator: 



Ut") H*) "► G%(t" ~ *') = Kl e-* /s ° (t "" t ' ) ^ [ l + 6{t" - t>) - 7 _6>(f - t") 



(70) 



where the matrix elements V\j define the unitary transformation which diagonalizes the hopping matrix /£, while f® 
are the corresponding eigenvalues. 

The Wick contractions are most conveniently defined and computed using a diagrammatic technique, i.e. by 
applying the Feynman rules shown in Fig. [4] Just like in the standard quantum field theory, one can prove that 
the corrections to the normalization factor Z^> exactly cancel out with the contribution of disconnected diagrams, 
order-by-order in perturbation theory. 

From the zero-th order diagram shown in Fig. [3] we readily re-obtain the unperturbed conditional probability 

*?(k/|ki) 



P? ) (kf\k i )= -Gi ikf (-t)Gi tki (t)=vl n e i ^ ht V^ t Vl s e- i f"/ ht V s ^ = G k/kt (i) 



(71) 



13 



a 



t" 
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-fL^ioV 3l {t"-t')p n 



2h 



FIG. 4: Feynman rules for the effective field theory for charge propagation in the macromolecule. On the left panel, we show 
the hole's Feynman propagator, on the center panel the effective interaction Vi, and on the right panel the effective interaction 
V 2 . 






FIG. 5: The diagrams involved in the leading-order correction to the conditional probability. The first two diagrams from the 
left-hand-side are called "self-energy" -type. The third diagram is called "crossing" -type and the last is a "tad-pole" diagram. 
Analog diagrams exist for "self-energy" and "tad-pole" diagrams, in which the vibron propagators are coupled to the hole 
propagator on the left. 



The different types 4 of diagrams which contribute to P^kyllq) are shown in Fig. p] We note that the first two 
of such diagrams contain a "self-energy" -type correction to one of the propagators. The third diagram contains the 
interaction of forward and backward propagating holes and will be called the "crossing" -type diagram. Finally, the 
last diagram is a "tad-pole" . After collecting all terms, we obtain the following expression for the first-oder correction 
to the conditional probability: 



^ (1) (kf|ki) 



4M7 

f3h 2 

2 

h 
2 
h 

2M7 
ph 2 



Re 



Im 
Im 



drdr' G k , kf (-*) G kfq ,(i - r') / q , s , A^r' - r) G° s , s (r' - r) / s ' q G^r) 



I drdr' G° kikt (-t) G° fq ,(t - r') / q , s , V 3i (r> - r) G% s {t> - r) /*, G qkl (r) 

f drdr' G kikt (-t) G° fq (i - r) / qs G s ° ki (r) Vy(r' - r) f s , s , 
Jo 

f drdr' G klS ,(t - r') /f, q , G q , kf (r') A j{ (r' - r) G kfS (t - r) /* q G qk ,(r) . 



(72) 



The first two lines are the contributions due to the "self-energy" -type diagrams, the third is derived from the "tad- 
pole" diagram, and the last line is derived from the "crossing" -type diagram. Further the details of this calculation 
are provided in Appendix [Aj 

We conclude this section by discussing the regimes in which we expect the perturbative expansion to be applicable. 
To this end, we introduce the explicit expressions for G^, A^-, and Vij into Eq. (72), take the short time limit 
t <C 1/fifc? t <C I/7, and consider for simplicity only the over-damped and under-damped asymptotic expressions for 
the Green's functions. 



Clearly, in addition to the diagrams shown in Fig. pi there are also equivalent ones in which the vibronic propagators are emitted and 
absorbed by the backward propagating fields. 
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FIG. 6: Upper panel: three-dimensional structure of a P3HT polymer. Lower panel: the coarse-grained representation corre- 
sponding to our effective model. 



In the in the over-damped regime, we find the conditions of validity 

<C 1 (from the Vi-type interaction) , 



2 m t 2 



p m n 2 k n 2 

4 m t 2 



< i 



(from the V-j-type interaction) . 



1)1 
In the under-damped regime, the conditions of validity of the perturbative expansion are 



2 4 S 2 e 

/? m n 2 k n 2 
2 As 2 1 2 

M flu h 



< 1 



< 1 



(from the Vi-type interaction) , 
(from the V-j-type interaction) . 



(73) 
(74) 

(75) 
(76) 



We note that in both the under-damped and over-damped regimes the perturbative approach is only valid at short 
times. This is completely expected, because the long-time and long-distance propagation necessarily involves multiple 
scattering of the hole with the heat bath and molecular vibrations, hence require a non- perturbative treatment. By 
plugging order of magnitude estimates of the normal- mode frequencies (fl*. ~ 10 -3 fs~ ), the gradient of the hopping 
matrix elements (/q ~ 10 -2 eV A -1 ), and the viscosity (7 = 0.1 fs~ ) we find that, at room-temperature and in the 
over-damped limit, the Vi-type interaction is several orders of magnitude larger than the V2-type, hence determines 
the range of validity of the perturbative expansion. On the other hand, in the opposite under-damped limit, the 
driving interaction is Va-type. 



V. CHARGE PROPAGATION THROUGH A POLY(3-ALKYLTHIOPHENE) MOLECULE 

Let us now illustrate the formalism developed in the previous sections by investigating the intra-chain propagation 
of electron holes through the backbone of a poly(3-alkylthiophene) (P3HT) polymer. Quasi-cristalline materials made 
of inter-digited PH 3T po lymers have received much attention, in connection with the possibility of realizing nano-scale 
organic transistor j 8 l 29 l 3Q [ The atomistic three-dimensional structure of a PH3T molecule is shown in the upper panel 
of Fig. |6| 

Here, we are only interested in providing a qualitative description of the charge propagation in such systems, leaving 
a more sophisticated quantitative description to our future work. Our main goal is to estimate the order-of-magnitude 
of the range of times and distances over which the perturbative approach is applicable. Secondly, we are interested in 
comparing the probability densities obtained by means of the perturbative calculation and by brute-force integration 
of the equation of motion (11). Finally, we investigate how the different charge-charge effective interactions which 



appear in Eq. (44) affect the charge dynamics and in particular quantum de-coherence and re-coherence phenomena. 
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FIG. 7: Left panel: DTP-B calculation of the molecular potential energy of two consecutive monomers as a function of the 
difference of dihedral angles. The mechanical equilibrium is retained at an angle s z which in this case is c± 20°. In general 
s l — +(— )0q for i odd (even), with Oq — 20°. Right panel: DFT-TB calculation of the transfer integrals Ti+i^ as a function of 
the difference in the dihedral angles (0i+i — Qi), in the proximity of the mechanical equilibrium configuretion 0q. 



Coarse- Grained Model 



In order to address these points it is sufficient to adopt a simple coarse-grained representation of the molecule, 
in which side-chain degrees of freedom are not taken explicitly into account. Furthermore, the molecular potential 
energy function is assumed to effectively depend only of the dihedral angles formed by neighboring aromatic rings. 
Hence, the molecular configuration is specified by the set of dihedral angles 9 = (#i, . . . , On) and the chain is mapped 
into an effective one-dimensional system consisting of N plaquettes which can rotate around their symmetry axis, as 
sketched in the lower-right panel of Fig. [6j 

The potential energy of a molecular configuration is approximated with sum of pairwise terms, each of which 
depends on the relative dihedral angle of two consecutive monomers, 



U(Q) = J2< 6 i~ 6 i+i)- 



(77) 



We have obtained the pair- wise interaction energy u(0i — #i+i) as a function of the relative angle = 6i — i+ i from 
DFT-TB electronic structure calculations, using the DFTB+ package 31 . The results are shown in the left panel of 
Fig. [7J In the lowest-energy configuration, the aromatic rings in the different residues form a relative dihedral angle 
of (—1)* #o (where 9$ = 20° and i is the monomer index). 

The conformational dynamics is therefore described by the Hamiltonian: 



H 



M 



1 N 

- y 

21 ^ 



P? + U(S) , 



(78) 



where I is the momentum of inertia of the monomers (including the contribution from the atoms in the side-chain) 
and pi = 10 i is the canonical momentum conjugated to the 0i generalized coordinate. 

At room temperature, this system performs only small thermal oscillations around the minimum-energy configura- 
tion. Hence, we can perform a small-angle expansion, leading to the simple harmonic form: 



N N-l 

h m ~ J- £ri + J2 ?(«<+! -<>i + (-i)^o) 2 + \e\ + f t 



(79) 



i=l 



2=1 



The last two terms follow from assuming that the first and last monomers in the chain are bond to external non- 
conducting leads which tend to align them horizontally. 

The momentum of inertia / can be calculated directly from the three-dimensional structure of the chain and affects 
the frequencies of the chain's normal modes of oscillations, 



(80) 
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TABLE I: Parameters or the coarse-grained model which describes intra-chain hole propagation in P3HT polymers. The PH3T 
chain investigated in the simulations consists of 9 monomers. 



e [eV] 


To [eV] 


Ti [eV] 


0o [deg] 


7 \fs~ L \ 


k [eV] 


K e ff [eV] 


T [°K] 


I [uma A 2 ] 


5.4 


0.4 


0.06 


20 


0.1 


0.13 


0.20 


300 


3400 



In the systems which are of technological and experimental interest, P3HT polymers are embedded in organic frame- 
works. In such a configuration, the chain exchanges energy with neighboring molecules, which play the role of a heat 
bath. In addition, the steric interaction with neighbors generates strong constraints on the chain dynamics and in 
particular affects the amplitude and frequencies of thermal oscillations. In order to account for this effect, we consider 
an effective model in which the spring constant k which appears in the molecular potential energy function U{&) is 
artificially rescaled in such a way that the typical square fluctuations of the dihedral angles around their equilibrium 
values, (A0 2 )md — jf X^(((^+i — @i + ( — 1)* 9o) 2 )md matches the value obtained from classical molecular dynamics 
simulations for a system of inter-digited PH3T polymers^ 



^eff ■ 



and 



(A6 2 ) 



MD 



k B T 

Keff 



(81) 



Also the hopping matrix elements T^+i between neighboring monomers and the on-site energies ei as a function of 
the chain's configurations have been obtained by DFT-TB calculations. In analogy with molecular energy calculations, 
the transition matrix elements T^+i have been computed assuming that they effectively depend only on the relative 
angles, \0 i+ i — 6i\. For sake of simplicity, we have taken the on-site energies e$ to be constant and equal to the value 
at the mechanical equilibrium configuration. This choice can be motivated by the observation made by different 
groups (see e.g. RefpJ) that fluctuations of the on-site energies have a much smaller effect on the electric conduction 
than fluctuations of the transfer matrix elements. The results for T^+i(^+i — #*), in the vicinity of the equilibrium 
configuration O are reported in the right panel of Fig. [7] By assuming a linear approximation, Eq. ([5| takes the form 

fmnifii) = f mn + f mn (\0 m — n \ - Oq) , (82) 



where 



fmn — ^0 (1 — $mn) ~ €-0 ^n 
fmn = ^1 (1 - 5mn) • 



(83) 
(84) 



Finally, the viscosity parameter 7 may be determined from MD simulations by computing the velocity autocorrelation 
function. On the other hand, we have observed that the results of the perturbative calculation depend very weakly 
on the this parameter. Hence, for sake of simplicity, here we assume a reasonable value 7 = 0.1 fs _1 . 

The numerical values of the parameters of this coarse grained model are summarized in Table [l| The equilibrium 
configuration can be chosen to be: 



f if 
I 0o if 



i odd, 
i even. 



(85) 



The hole propagator G® m (t) is constructed by diagonalizing the f° matrix, defined in Eq. (83). Its n-th eigenvector 
reads 



<t>n{j) 



A/" + l 



sin 



The corresponding eigenvalue is 



where fc„ is the wave- vector 



(N + l) 
E„ = -e - 2T cos(k n ) , 



(86) 
(87) 



k n 



Tin 



n = 1,2,. ..,iY. 



(N+l)' 
Hence, the hole Feynman propagator is given by 

2(N+1) 

G UW = E PnVf) Mi) e-i^ ( 7 + 6>(*)-7-0(-*)). (89) 

n=l 

With the present set of model parameters, we find that the friction coefficient 7 is significantly larger than the typical 
frequencies of normal models, i.e. 7 ^> Qk- Hence, we can use the over-damped limit for the vibronic propagators. 
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FIG. 8: Time evolution of the charge density at the right endpoint of the chain, assuming that a hole is initially created at 
left-end of the chain. The solid line represents the unperturbed prediction, the dashed line includes the effects of the coupling 
with the molecular vibration and the heat bath to leading-order in perturbation theory. 



B. Time Evolution of the Conditional Probability 

We assume that the hole is initially created at the left end of the chain and evolves in time at a temperature of 
300 K. In Fig. [8] we show the results of our leading-order perturbative calculation of the hole density at the opposite 
end of the chain as a function of the time interval £, i.e. Pt(kzv|ki), where N = 9. 

Some comment on these results are in order. First of all, we notice the existence of three peaks in the charge density, 
at times t ~ 10, 20, 40 ps. These corresponds to integer multiples of the time interval it takes the hole to run along 
the entire chain. Next, we note that the scattering of the holes with the molecular normal modes and with the heat 
bath slows down the charge propagation, as expected. This is evident from the fact that the probability of observing 
the hole at the right end-point of the chain as a function of time is reduced once the perturbative corrections are 
included. Correspondingly, the times at which the hole rebounds to the right end-point of the chain are delayed by 
the interactions. We recall that the norm of the conditional probability is conserved up to corrections which are of 
higher-order in the perturbative expansion. Hence, the reduction of the charge density at the end-point of the chain 
implies that the scattering distributes the charge density in the central region of the chain. 

We observe that the correction to the conditional probability starts to be of the same order of the unperturbed 
prediction starting from time intervals of the order of ~ 40 fs. Beyond this time scale, the perturbative approach 
breaks down and one has to resort on non-perturbative approaches in which many Feynman diagrams are re-summed. 



C. Quantifying the Loss of Quantum Coherence 



We have seen that analytic perturbative calculations break down beyond a few tens of fs, hence do not represent a 
useful tool to investigate the long-time long-distance dynamics of hole propagation. On the other hand, they provide a 
valuable tool to gain analytic insight into the physical mechanisms which drive de-coherence and re-coherence during 
hole propagation across the chain. „„ 

As measure of the degree of quantum coherence in the dynamics of an open system, we consider the raticP^™ 



R(t) 



Tr p 2 (t) 
Trp(t) 



(90) 



In Appendix [B] we show that this ratio is identically equal 1 for pure states (corresponding to fully coherent propa- 
gation), and that it is smaller than 1 for mixed states. 

In Fig. [9] we compare this ratio for the model under consideration in the limit of unperturbed propagation and 
including the leading-order scattering with the molecular vibrations and with the heat bath. We see that the interaction 
with the environment suppresses the quantum coherence on time scales which are of the order of 10 fs. 

It is interesting to compare the contribution to R(t) coming from the different Feynman diagrams shown in Fig. [51 
We find that the quantum decoherence is driven by the "cross" -type diagram shown in figure [5j which tends to 
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FIG. 9: Perturbative calculation of the different contribution to the ratio R(t) which quantifies the de-coherence effects in 
hole propagation. The dot-dashed line is the unperturbed result, which remains coherent at all times. The solid line is the 
result including the leading order correction. Dotted and dashed lines correspond to the contribution of "self-energy" -type and 
"crossing" -type, respectively. For this model, the "tad-pole" -type contribution is identically zero. 



correlate the field components associated to propagation forward and backwards in time. In the equivalent zero- 
temperature quantum-field theory effective picture, the quantum decoherence emerges as a result of the formation 
of a particle- antiparticle "bound state". On the other hand, the so-called "self -type diagrams act in the opposite 
direction, slowing down the overall rate of quantum decoherence. 

The identification on the diagram which drives the quenching of R(t) with time offers a scheme to study how 
the chemical and mechanical properties of the macromolecule affect the quantum decoherence of the propagating 
excitation. Indeed, by varying the parameters of the effective theory (namely the spectrum of normal modes Uk 
entering in the vibron propagators, and unperturbed tight-binding matrix elements f® m and their gradient /j* m , which 
enter the effective interaction vertexes) and computing the corresponding relative weight of the "cross" -type diagram, 
one may in principle identify what properties of the macromolecular system are most effective in suppressing (or 
enhancing) the quantum coherent transport. This information may be useful e.g. in the context of the study of 
exciton propagation in photosynthetic complexes, which have been found to display quantum coherent dynamics over 
surprisingly long time intervals. 

D. Comparison between the perturbative estimate and the result of integrating the quantum/stochastic 

equations of motion 



The perturbative approach developed in the previous sections allows to analytically compute the charge density, 
in the range of time intervals < t < 50 fs. It is interesting to compare the perturbative calculation with the 
results of non-pert urbative numerical simulations, obtaine d by averaging over many independent solutions of the set 
of quantum/stochastic equations of motion defined in Eq. ( pT| ) and derived in RefP^. On the one hand, this provides 
a non-trivial test for the perturbative scheme developed in this work. On the other hand, it offers an estimate of 
the statistical accuracy which is needed in order to resolve the effects of the interactions on the charge propagation 
dynamics in non-perturbative numerical simulations. 

In Fig. [To] we present the difference between the interacting and the free conditional probabilities P t (N,l) — 
Pt (N, 1), evaluated in the perturbative and non-perturbative methods (we recall that N = 9). The shaded area 
represents the statistical error in numerical simulations, which is estimated from the variance calculated from 10000 
independent trajectories. Accumulating this statistics required about 6 Central Processor Units (CPU) hours of 
simulation on a regular desktop. By contrast, the perturbative estimates took about a minute on the same machine. 

We find that the two approaches are quantitatively consistent with one another, even at time scales of the order of 
50 fs. Beyond such a time scale the perturbative approach becomes unreliable and the comparison is meaningless. 
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FIG. 10: Time evolution of the charge probability density at the right end-point of the chain assuming that the hole is initi ally 
created at left-end of the chain. We compare non-perturbative numerical calculation obtained by integrating the Eq.s (11) 
(solid line) and the analytic perturbative calculations (dotted line). The shaded area represents the statistical uncertainty on 
the non-perturbative calculation. 



It is important to emphasize that these two methods are based on different approximations. In particular, the 
algorithm defined by Eq.s (11) was obtained by neglecting the fluctuations of the coherent fields around their func- 



tional saddle-point solution. At such a saddle-point, the forward- and backward- propagating fields are identical, 
(j)'{t) = (j)"{t). The leading-order perturbative estimate goes beyond such a saddle-point condition and accounts for 
independent quadratic fluctuations on 4> r (t) and on <j)" (t). The relatively good agreement between the two calculation 
schemes at short times can be used as an argument in favor of the accuracy of the saddle-point approximation used 
in the non-perturbative approach. 



VI. CONCLUSION AND OUTLOOK 



In this work we have introduced a path integral approach to investigate the real-time dynamics of a quantum 
excitation propagating through a macromolecular system. In our starting model, all the atomic coordinates are 
explicitly taken into account, the dynamics of the quantum excitation propagating throughout the molecule is described 
by an effective configuration-dependent tight-binding Hamiltonian, and the heat bath is represented by the Caldeira- 
Leggett model. 

By adopting the first quantization for the molecular and the heat bath variables, and the second-quantization for the 
quantum excitation variables, we have been able to analytically trace out from the density matrix both the heat bath 
and the atomic nuclear dynamics. The result is a path integral formulated solely in terms of the quantum excitation's 
degrees of freedom. In contrast to other field-theoretic approaches to charge transport in macromolecules 3 ^^, in the 
present effective theory approach fluctuation-dissipation effects are fully taken into account. In addition, the free 
parameters can be obtained directly from microscopic quantum chemistry calculations. 

A first important advantage of our effective field theory formalism is that it makes it possible to describe real- 
time dynamics in an open quantum system without employing the Keldysh contour: indeed, the degrees of freedom 
associated to the backwards time-propagation in the Keldysh formalism can be effectively re-interpreted as anti-matter 
components of forward propagating quantum fields. In this way, the path integral which yields the reduced density 
matrix for the quantum excitation becomes formally equivalent to one associated to a standard vacuum-to-vacuum 
two-point correlation function, in a closed system at zero-temperature. This formal analogy is quite useful, as it 
immediately yields the Feynman rules to compute by perturbation theory the corrections to the density matrix due 
to the interaction between the propagating excitation, the atomic coordinates and the heat bath degrees of freedom. 

The possibility of performing analytic calculations in the weak-coupling and short-time regime opens the door to 
a detailed investigation of the effects which generate quantum decoherence in molecular systems coupled to a heat 
bath. We have identified a specific Feynman diagram which dominates the dissipation of the quantum coherence, 
and correlates forward and backward propagating fields. In the future, it would be interesting to perform the same 
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calculation for a class of simple models, in order to clarify whether the dominance of this diagram is a universal 
property of all open quantum systems, and if so unveil its physical interpretation. 

For illustration purposes, we developed a coarse-grained model and applied this formalism to investigating the 
intra-chain propagation of holes in a P3HT polymer. We found that the propagation can be described in perturbation 
theory up to about 40-50 fs. 

Beyond that time scale, non-perturbative approaches are required. Since the coherent-state path integral is af- 
fected by a dynamical sign problem, Monte Carlo approaches would be challenging. An alternative non-perturbative 
approach consists in directly integrating the quantum and stochastic equation of motions (11) which follow from a 



functional saddle-point approximation. The comparison with the analytic perturbative calculations has shown that 
the underlying saddle-point approximation is quite robust. However, a large statistics seems to be necessary in order 
to resolve the tiny effects of the interaction with the heat bath and with the vibronic modes. An alternative non- 
perturbative approach which would not involve stochastic averages and MD simulations could be provided by the 
self-consistent saddle-point approximation of Eq. (44). We plan to develop and test such an approach in our future 
work. 

Finally, an obvious direction for our future research consists in developing the formalism to compute the electric 
current induced by an external electric field. 
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Appendix A: Details on the vibronic Green's functions structure and the perturbative calculations 



In performing the path integral over the Sr and y variables, we exploit the standard result for Gaussian functional 
integrals: 



fvcj) exp - / dt'dt" ^it^Aijit' -t")(j) j {t")+ f dt' Bi(t')^i(t') + C 
J I Jo Jo 

\ [ dt'dt" Bi(t') At 1 it' - t") B 3 (t") + C 
4 Jo 



detA- 



(Al) 



The vibronic two-point functions Aij(t f — t") and Vij(t — £'), which enter in Eq. (52) are contracted from the Green's 
functions of the L^L, L and L^ operators: 



Ay(t) 



DL] (t) = [ (MdfSij + jMdtSij + Hij) (MdfSij - jMdtSij + Hij) ] * , 



(A2) 



D 



(*)■ 



L 



(t) = [M (ft 2 + 7 ft) Sij + Hij] X + [M (d 2 t - 7 ft) 5tj + Hij] 1 . (A3) 



In order to compute them it is convenient to consider the Fourier transform to frequency space. We also transform 
into the normal mode basis, by applying the unitary transformation U which diagonalizes the Hessian operators. We 
obtain: 



M = T72 U in [ (^ 2 - ^ - fin) (^ 2 + ilU ~ Q n ) ] 



M 2 



U n 



(uj 2 — ijuj — ft n ) + (uj 2 + ijuj — fi n ) 



U nj •> 



(A4) 
(A5) 



where Hij is the Hessian, and Q n are the corresponding normal modes. The expressions (54) and (55) for A^-(t) 
and Vij(t) are obtained by Fourier transforming back to the time representation, taking the continuumTimit for the 
Fourier sum. 
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The following traces enter the derivation of the perturbative estimate (72): 



tr [7-7V7 5 ] = N spin (A6) 

-tr [7-7 ] = tr [7+7 ] = N spin (A7) 

-tr [7-7V7V7 5 ] = tr [7"7 5 7 + 7°7 + 7 5 ] = N spin (A8) 

tr [tVtVtVtV] = tr [ 7+7 5 7 - 7 o 7 - 7 o 7 - 7 5J = ^^ (A9) 

where N spin is the degeneracy number associated to the spin of the excitation and is equal to 2 for spin- 1/2 fermions 
and 1 for spin-0 bosons. We recall that also the dimensionality and definition of the 7-matrixes depend on the spin 
of the propagating particle. In particular, for spin-0 bosons they reduce to 70 = 73, and 75 = ti, where T3 and t\ are 
Pauli matrixes, 

Appendix B: Measuring Quantum Decoherence 

In this appendix we review the proof that the ratio 

R(t)=Tr[p 2 (t)}/Tr[p(t)}, (Bl) 

provides a measurement of the degree of decoherence of the system. 

We first consider a pure state, denoted by the vector |x), and we represent the corresponding density operator with 
P=\X)(X\, sothatTV[ / 5] = ( X |x). 

The operator p 1 reads p 2 = |x)(xlx)(xl> while its trace is 

Tr[p 2 ] = (xlx) = (Tr[p]) 2 , (B2) 

hence R(t) = 1. For a mixed state, there is no single state vector describing the system, and R(t) < 1. 
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